Figure 1: code for insets

Coverage example

## Project population
vacc_reconstruct <- function(t, y, parms){
  dS <- -parms["deaths"] * y["S"] + parms["births"] * (y["S"] + y["V"]) + parms["waning"] * y["V"]
  dV <- -parms["deaths"] * y["V"] - parms["waning"] * y["V"]
  dWaning <- -parms["waning"] * y["V"]
  dDeaths <- -parms["deaths"] * y["V"]
  dBirths <- -parms["births"] * (y["S"] + y["V"])
  return(list(c(dS, dV)))
} 

eventfun <- function(t, y, parms){
  with (as.list(y),{
    diff <- 0.7*(S + V) - V ## newly vaccinated
    V <- V + diff
    S <- S - diff
    return(c(S, V))
  })
}


## times in relevant timestep
times <- seq (1/52, 5, by = 1/52) ## weekly for 10 years

## initial pop
IS <- 70000
IV <- 0

## param
births <- 0.50
deaths <- 0.42
waning <- 0.33
parms <- c(deaths = deaths, births = births, waning = waning)

## run
calcsus <- lsoda(y = c(S = IS, V = IV), times = times, func = vacc_reconstruct, parms = parms, 
                 events = list(func = eventfun, time = seq(1/52, 5, by = 1)))
vacc_df <- data.frame(Year = calcsus[, "time"], propVacc = calcsus[,"V"]/(calcsus[,"V"] + calcsus[,"S"]))
vacc_a <- ggplot(data = vacc_df, aes(x = Year, y = propVacc)) +
  geom_line(color = "lightseagreen", size = 1) +
  xlab ("Year") +
  ylab("Coverage") + 
  ylim(c(0, 1)) +
  geom_hline(yintercept = 0.69, color = "aquamarine4", linetype = 2, size = 1) +
  newtheme +
  theme(text = element_text(color = "aquamarine4"), axis.text = element_text(color = "aquamarine4"), 
        axis.line = element_line(color = "aquamarine4"))
ggsave("figs/fig1/vacc_a.jpeg", vacc_a, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/vacc_a.jpeg")

Secondary cases

## params
secondaries <- data.frame(secondaries = rnbinom(1000, mu = 1.2, size = 1.3))
sec_b <- ggplot(data = secondaries, aes(x = secondaries)) +
  geom_histogram(binwidth = 1, color = "lightgray", fill = "red", alpha = 0.5) +
  geom_vline(xintercept = 1.2, color = "red", linetype = 2, size = 1.1) +
  xlab("Secondary cases") +
  ylab("Frequency") +
  newtheme +
  theme(text = element_text(color = alpha("firebrick3", 1)),
        axis.text = element_text(color = alpha("firebrick3", 1)), 
        axis.line = element_line(color = alpha("firebrick3", 1)))
ggsave("figs/fig1/sec_b.jpeg", sec_b, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/secc_b.jpeg")

Incubation period

incubation <- as.data.frame(list(Days = seq(0, 365, 1), 
                           Density = dgamma(seq(0, 365, 1), shape = 1.1518, 
                                    rate = 0.0429)))
inc_c <- ggplot(data = incubation, aes (x = Days, y = Density)) +
  geom_line(color = "darkred", size = 1.2) +
  xlab("Incubation Period \n (days)") +
  geom_vline(xintercept = 22.3, color = "darkred", linetype = 2, size = 1.1, alpha = 0.5) +
  newtheme +
  theme(text = element_text(color = "darkred"), axis.text = element_text(color = "darkred"), 
        axis.line = element_line(color = "darkred"))
ggsave("figs/fig1/inc_c.jpeg", inc_c, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/inc_c.jpeg")

Infectious period

infectious <- as.data.frame(list(Days = seq(0, 10, 0.1), 
                           Density = dgamma(seq(0, 10, 0.1), shape = 2.9012, 
                                    rate = 1.013)))
inf_d <- ggplot(data = infectious, aes (x = Days, y = Density)) +
  geom_line(color = "navy", size = 1.2) +
  xlab("Infectious Period \n (days)") +
  geom_vline(xintercept = 3.1, color = "navy", linetype = 2, size = 1.2, alpha = 0.5) +
  newtheme +
  theme(text = element_text(color = "navy"), axis.text = element_text(color = "navy"), 
        axis.line = element_line(color = "navy"))
ggsave("figs/fig1/inf_d.jpeg", inf_d, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/inf_d.jpeg")

Road network with scale

tz_roads <- readOGR("data/Tanzania_Roads/Tanzania_Roads.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/mrajeev/Documents/Projects/ModelingChapter/data/Tanzania_Roads/Tanzania_Roads.shp", layer: "Tanzania_Roads"
## with 518 features
## It has 17 fields
tz_roads <- fortify(tz_roads)

map_e1 <- ggplot(data = tz_roads, aes(x = long, y = lat, group = group)) +
  geom_line(color = "mediumorchid4") +
  scalebar(tz_roads, location = "bottomleft", 
           dist = 200, dist_unit = "km", transform = TRUE, box.fill = c("mediumorchid4", "grey"),
           box.color =  c("mediumorchid4", "grey"),
           st.size = 5, st.dist = 0.05, st.color = "mediumorchid4", 
           height = 0.02, model = 'WGS84') +
  theme_map()

ggsave("figs/fig1/map_e1.jpeg", map_e1, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/map_e1.jpeg")

Dispersal kernel

## params
dispersal <- as.data.frame(list(km = seq(0, 10, 0.01), 
                           Density = dgamma(seq(0, 10, 0.01), shape = 0.215, 
                                    scale = 0.245)))
disp_e2 <- ggplot(data = dispersal, aes (x = km, y = Density)) +
  geom_line(color = "mediumorchid4", size = 1.2) +
  xlab("Distance to next bite \n (km)") +
  geom_vline(xintercept = 0.88, color = "mediumorchid4", linetype = 2, size = 1.2, alpha = 0.75) +
  newtheme +
  theme(text = element_text(color = "mediumorchid4"), axis.text = element_text(color = "mediumorchid4"), 
        axis.line = element_line(color = "mediumorchid4"))

ggsave("figs/fig1/disp_e2.jpeg", disp_e2, device = "jpeg", width = 5, height = 5)
# include_graphics("figs/fig1/disp_e2.jpeg")

Final figure 1

include_graphics("figs/fig1.jpeg")

Figure 2: Summarizing studies

studies <- read.csv("data/modeling_studies.csv")

## A. Studies by Year
year_A <- ggplot(data = studies, aes(x = Year)) +
  geom_bar(fill = "grey50") +
  scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
  ylab("Number of studies") +
  labs(tag = "A") +
  newtheme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## B. Studies by country
## refactor countries first
studies %>%
  mutate(Country = fct_recode(Country, Other = "India", Other = "Israel", Other = "Kenya")) -> studies
## then plot   
country_B <- ggplot(data = studies, aes(x = fct_relevel(reorder(Country, Country, function(x)-length(x)), 
                                            "Multiple (Africa)", "Multiple (global)", 
                                           "Other", "Hypothetical", after = 10))) + 
  geom_bar(fill = "grey50") +
  scale_y_continuous(breaks = seq(0, 10, by = 2), labels = scales::number_format(accuracy = 1)) +
  ylab("Number of studies") +
  xlab("") +
  labs(tag = "B") +
  newtheme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## C. R0 plot
r_ests <- read.csv("data/r_ests.csv")
R0_C <- ggplot(data = filter(r_ests, R_type == "R0"), aes(x = Estimate)) +
  geom_histogram(binwidth = 0.25, color = "white", fill = "grey50") +
  scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
  labs(tag = "C") +
  xlab("R0 estimate") +
  ylab("Number of studies") +
  newtheme

## D. Key modeling decisions and features
studies %>%
  select(14:20) %>%
  gather() %>%
  filter(value %in% "Y") %>%
  group_by(key) %>%
  summarize(proportion = n()/51) -> summary
features_D <- ggplot(data = summary, aes(x = reorder(key, -proportion), y = proportion)) +
  geom_col(fill = "grey50") +
  scale_x_discrete(labels = c("DDT", 
                              "Fit to data", "Stochastic", "Spatial", 
                              "Heterogeneity", "Introductions", 
                              "Underreporting")) +
  ylab("Proportion of studies") +
  ylim(c(0, 0.8)) +
  xlab("") +
  labs(tag = "D") +
  newtheme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))


fig2 <- year_A + R0_C + country_B + features_D + plot_layout(ncol = 2, nrow = 2)
ggsave("figs/fig2.jpeg", fig2, device = "jpeg", width = 20, height = 15)
## Warning: Removed 1 rows containing non-finite values (stat_count).
include_graphics("figs/fig2.jpeg")

Figure 3: data

## Reported data
repdata <- read.csv("data/reported_data.csv")
repdata %>%
  mutate(Nobs_maxed = case_when(Nobs < 250 ~ as.numeric(Nobs), Nobs >= 250 ~ 250), 
         Spatial.scale = fct_relevel(Spatial.scale, "National", after = Inf), 
         Temporal.scale = fct_relevel(Temporal.scale, "Month", after = 2)) -> repdata

data_A <- ggplot(data = repdata, aes(x = reorder(Type.of.data, Type.of.data,function(x)-length(x)))) + 
  geom_bar() +
  newtheme +
  scale_x_discrete(labels = c("Clinical and lab \n confirmed animal cases", 
                              "Reported human \n deaths", 
                              "Sequence data", "Lab confirmed \n animal cases", 
                              "Lab confirmed \n human exposures", "Reported \n animal bites")) +
  xlab("Type of data") +
  ylab("Number of data \n sources") +
  labs(tag = "A") +
  coord_flip()

data_B <- ggplot(data = repdata, aes(x = Temporal.scale, y = Length, fill = Spatial.scale, size = Nobs_maxed)) +
  ggbeeswarm::geom_beeswarm(cex = 3.5, alpha = 0.75, shape = 21, color = "grey50", stroke = 1.5) +
  scale_fill_manual(values = c("Darkred", "Red", "Pink", "Grey50"), name = "Spatial scale") +
  scale_size_area(name = "# of Observations", breaks = c(10, 100, 150, 250), 
                       labels = c("10", "100", "150", "250+")) +
  ylab("Length of time \n series (years)") +
  xlab("Temporal scale") +
  newtheme +
  labs(tag = "B") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
  guides(fill = guide_legend(override.aes = list(size = 3)))
fig3 <- data_A + data_B + plot_layout(ncol = 2, widths = c(1, 1))
ggsave("figs/fig3.jpeg", plot = fig3, device = "jpeg", height = 6, width = 15)
include_graphics("figs/fig3.jpeg")

Supplementary Figure S1

# models
SEI.rabies <- function(nu = 0.50, mu = 0.42, R0 = 1.2, sigma = (1/22.3*365),
                       gamma = (1/3.1*365), K = 20, START.S = 50000, 
                       START.E = 0, START.I = 2,
                       START.S.density = 15, START.I.density = 0.01, 
                       START.E.density = 0.00, years = 20, steps = 1/52, model = "frequency",...){
  
  ### Deterministic skeleton of models for rabies
  
  require(deSolve)
  
  times <- seq(from = 0, to = years, by = steps)
  
  ## Simple FDT model with complete disease-induced mortality
  if (model == "frequency"){
    
    START.N <- START.S + START.E + START.I
    
    beta = (R0*gamma*(sigma + mu))/sigma
      
      # This function models a time step for the SIR:
      dx.dt.SIR <- function(t, y, parms) {
        N <- y["S"] + y["E"] + y["I"] 
        S <- y["S"]
        E <- y["E"]
        I <- y["I"]
        
        # Calculate the change in Susceptible
        dS <-  parms["nu"]*(S + E) - parms["mu"]*S - parms["beta"]*S*I/N
        
        #Calculate the change in Exposed
        dE <- parms["beta"]*S*I/N - parms["mu"]*E - parms["sigma"]*E
        
        # Calculate the change in Infected
        dI <- parms["sigma"]*E - parms["gamma"]*I
        
        # Return a list with the changes in S, E, I, R at the current time step
        return(list(c(dS, dE, dI)))
      }
      
      # Create the parameter vector
      parms <- c(nu = nu, mu = mu, beta = beta, sigma = sigma, gamma = gamma)
      inits <- c(S = START.S, E = START.E, I = START.I)
  }
  
  ## Classical fox rabies model (widely used for dog rabies as well, adapted from Anderson et al. 1981)
  if (model == "density"){
    
    START.N.density <- START.S.density + START.I.density + START.E.density
    
    # the old school rabies model
    # This function models a time step for the SIR:
    dx.dt.SIR <- function(t, y, parms) {
      N <- y["S"] + y["I"] + y["E"]
      S <- y["S"]
      E <- y["E"]
      I <- y["I"]

      ## Calculate dmort and beta here from parms and plug them in
      dmort <- (parms["nu"] - parms["mu"])/parms["K"]
      beta <- (parms["R0"]*(parms["sigma"] + parms["nu"])*(parms["gamma"] + parms ["nu"]))/(parms["sigma"]*parms["K"])
      
      # Calculate the change in Susceptible
      dS <-  parms["nu"]*N - parms["mu"]*S - dmort*S - beta*S*I
      
      #Calculate the change in Exposed
      dE <- beta*S*I - parms["mu"]*E - dmort*N*E - parms["sigma"]*E
      
      # Calculate the change in Infected
      dI <- parms["sigma"]*E - parms["mu"]*I - dmort*N*I - parms["gamma"]*I
      
      # Return a list with the changes in S, E, I, R at the current time step
      return(list(c(dS, dE, dI)))
    }
    
    # Create the parameter vector
    parms <- c(nu = nu, mu = mu, R0 = R0, K = K, sigma = sigma, gamma = gamma)
    inits <- c(S = START.S.density, E = START.E.density, I = START.I.density)
  }
  
  # Run the ODE solver
  SIR.output <- lsoda(y = inits, 
                      times = times, 
                      func = dx.dt.SIR, 
                      parms = parms)
  SIR.output <- as.data.frame(SIR.output)
  SIR.output$N <- SIR.output$S + SIR.output$E + SIR.output$I
  return (SIR.output)
}

## Density dependent
ddt_1.01 <- SEI.rabies(R0 = 1.01, years = 20, model = "density")
ddt_1.01 <- data.frame(infected = unname(tapply(ddt_1.01$I, (seq_along(ddt_1.01$I)-1) %/% 4, sum)[1:260]),
                       pop = ddt_1.01$N[seq(1, 20*52, 4)], R0 = 1.01, trans = "Density")
ddt_1.1 <- SEI.rabies(R0 = 1.1, years = 20, model = "density")
ddt_1.1 <- data.frame(infected = unname(tapply(ddt_1.1$I, (seq_along(ddt_1.1$I)-1) %/% 4, sum)[1:260]),
                       pop = ddt_1.1$N[seq(1, 20*52, 4)], R0 = 1.1, trans = "Density")
ddt_1.05 <- SEI.rabies(R0 = 1.05, years = 20, model = "density")
ddt_1.05 <- data.frame(infected = unname(tapply(ddt_1.05$I, (seq_along(ddt_1.05$I)-1) %/% 4, sum)[1:260]),
                      pop = ddt_1.05$N[seq(1, 20*52, 4)], R0 = 1.05, trans = "Density")

ddt <- bind_rows(ddt_1.01, ddt_1.1, ddt_1.05)
ddt$time <- 1:260

ddt_plot <- ggplot(data = ddt, aes(x = time, y = infected/pop, group = R0, color = as.factor(R0))) + 
  geom_line(size = 1) +
  scale_color_manual(values = c("lightcoral", "red", "darkred"), name = "R0") +
  xlab("Months") +
  ylab("Incidence (monthly proportion \n of population infected)") +
  newtheme +
  labs(tag = "B", subtitle = "Density")

## Frequency dependent
fdt_1.01 <- SEI.rabies(R0 = 1.01, years = 20, model = "frequency")
fdt_1.01 <- data.frame(infected = unname(tapply(fdt_1.01$I, (seq_along(fdt_1.01$I)-1) %/% 4, sum)[1:260]),
                       pop = fdt_1.01$N[seq(1, 20*52, 4)], R0 = 1.01, trans = "Frequency")
fdt_1.1 <- SEI.rabies(R0 = 1.1, years = 20, model = "frequency")
fdt_1.1 <- data.frame(infected = unname(tapply(fdt_1.1$I, (seq_along(fdt_1.1$I)-1) %/% 4, sum)[1:260]),
                      pop = fdt_1.1$N[seq(1, 20*52, 4)], R0 = 1.1, trans = "Frequency")
fdt_1.05 <- SEI.rabies(R0 = 1.05, years = 20, model = "frequency")
fdt_1.05 <- data.frame(infected = unname(tapply(fdt_1.05$I, (seq_along(fdt_1.05$I)-1) %/% 4, sum)[1:260]),
                       pop = fdt_1.05$N[seq(1, 20*52, 4)], R0 = 1.05, trans = "Frequency")
fdt <- bind_rows(fdt_1.01, fdt_1.1, fdt_1.05)
fdt$time <- 1:260

fdt_plot <- ggplot(data = fdt, aes(x = time, y = infected/pop, group = R0, color = as.factor(R0))) + 
  geom_line(size = 1) +
  scale_color_manual(values = c("lightcoral", "red", "darkred"), name = "R0") +
  xlab("Months") +
  ylab("Incidence (monthly proportion \n of population infected)") +
  guides(colour = "none") +
  newtheme +
  labs(tag = "A", subtitle = "Frequency")
figS1 <- fdt_plot / ddt_plot
ggsave("figs/figS1.jpeg", plot = figS1, device = "jpeg", height = 11, width = 11)
include_graphics("figs/figS1.jpeg")

Supplementary Figure S2

## Study category
studies %>%
  mutate_at(vars(starts_with("Category")), 
            function(x) ifelse(x == "N" | is.na(x), 0, 1)) %>%
  select(Reference, starts_with("Category")) -> categories
names(categories)
## [1] "Reference"                             
## [2] "Category..explain.observed.patterns"   
## [3] "Category..estimate.key.parameters"     
## [4] "Category..study.control.measures"      
## [5] "Category..identify.drivers.of.dynamics"
## [6] "Category..predict.future.trends"
names(categories) <- c("Reference", "Explain patterns",
                       "Estimate parameters", "Study control",
                       "Identify drivers", "Predict trends")

## Upset plot
jpeg("figs/figS2.jpeg", width = 500, height = 500)
upset(categories, nsets = 5, nintersects = NA,
      sets.x.label = "Type of study", mainbar.y.label = "Number of studies", text.scale = 2)
dev.off()
## quartz_off_screen 
##                 2
## Bar plot
studies %>%
  select(Reference, starts_with("Category")) %>%
  gather() %>%
  filter(value %in% "Y") -> cat_summary
## Warning: attributes are not identical across measure variables;
## they will be dropped
table(cat_summary$key)
## 
##      Category..estimate.key.parameters 
##                                     15 
##    Category..explain.observed.patterns 
##                                     12 
## Category..identify.drivers.of.dynamics 
##                                     10 
##        Category..predict.future.trends 
##                                      1 
##       Category..study.control.measures 
##                                     26
include_graphics("figs/figS2.jpeg")

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] UpSetR_1.4.0         patchwork_1.0.0.9000 forcats_0.4.0       
##  [4] stringr_1.4.0        dplyr_0.8.3          purrr_0.3.2         
##  [7] readr_1.3.1          tidyr_1.0.0          tibble_2.1.3        
## [10] tidyverse_1.2.1      ggsn_0.5.0           ggthemes_4.2.0      
## [13] rgdal_1.4-6          sp_1.3-1             deSolve_1.24        
## [16] ggplot2_3.2.1        knitr_1.25          
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.1         jsonlite_1.6       modelr_0.1.5      
##  [4] assertthat_0.2.1   vipor_0.4.5        cellranger_1.1.0  
##  [7] yaml_2.2.0         pillar_1.4.2       backports_1.1.5   
## [10] lattice_0.20-38    glue_1.3.1         digest_0.6.22     
## [13] rvest_0.3.4        colorspace_1.4-1   htmltools_0.4.0   
## [16] plyr_1.8.4         pkgconfig_2.0.3    broom_0.5.2       
## [19] haven_2.1.1        scales_1.0.0       jpeg_0.1-8        
## [22] ggmap_3.0.0        generics_0.0.2     ellipsis_0.3.0    
## [25] withr_2.1.2        lazyeval_0.2.2     cli_1.1.0         
## [28] magrittr_1.5       crayon_1.3.4       readxl_1.3.1      
## [31] maptools_0.9-5     evaluate_0.14      nlme_3.1-141      
## [34] xml2_1.2.2         foreign_0.8-72     class_7.3-15      
## [37] beeswarm_0.2.3     tools_3.6.1        hms_0.5.1         
## [40] RgoogleMaps_1.4.4  lifecycle_0.1.0    munsell_0.5.0     
## [43] compiler_3.6.1     e1071_1.7-2        rlang_0.4.1       
## [46] classInt_0.4-1     units_0.6-4        rstudioapi_0.10   
## [49] rjson_0.2.20       bitops_1.0-6       labeling_0.3      
## [52] rmarkdown_1.16     gtable_0.3.0       DBI_1.0.0         
## [55] R6_2.4.0           gridExtra_2.3      lubridate_1.7.4   
## [58] zeallot_0.1.0      KernSmooth_2.23-15 stringi_1.4.3     
## [61] ggbeeswarm_0.6.0   Rcpp_1.0.2         vctrs_0.2.0       
## [64] sf_0.8-0           png_0.1-7          tidyselect_0.2.5  
## [67] xfun_0.10